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Abstract 

Based on experimental observations in Daphnia, we introduce an agent-based model for 
the motion of single and swarms of animals. Each agent is described by a stochastic equation 
that also considers the conditions for active biological motion. An environmental potential 
further reflects local conditions for Daphnia, such as attraction to light sources. This model 
is sufficient to describe the observed cycling behavior of single Daphnia. To simulate vortex 
swarming of many Daphnia, i.e. the collective rotation of the swarm in one direction, we ex- 
tend the model by considering avoidance of collisions. Two different ansatzes to model such 
a behavior are developed and compared. By means of computer simulations of a multi-agent 
system we show that local avoidance - as a special form of asymmetric repulsion between 
animals - leads to the emergence of a vortex swarm. The transition from uncorrelated rota- 
tion of single agents to the vortex swarming as a function of the swarm size is investigated. 
Eventually, some evidence of avoidance behavior in Daphnia is provided by comparing ex- 
perimental and simulation results for two animals. 
Keywords: active motion, swarming, zooplankton, Brownian agents 

1 Introduction 

Swarming is a prominent example of complex behavior in biological systems. This form of collec- 
tive motion may emerge from the interplay of individual behavior and local interactions of a large 
number of individuals (agents). Swarms (also called herds, flocks, schools) can often be observed 
in certain mammals, fish, insects, and birds for various benefits, such as enhanced feeding and 
mating as well as more successful predator avoidance [see e.g. SJ3 J2I El, IU Q El} . This has 



been reported for several prey animals, e.g. in planktivore fish [2l|, [46j, in some species of birds 
[3], as well as in zooplankton 2f], 3o| . 



Detailed experimental investigations on swarming, however, are rare, either because of the size 
of the animals or because well defined conditions for experiments are difficult to realize. Earlier 



chance observations were reported for horizontally circling zooplankton in the field 34j[. They 
triggered further experiments with Daphnia relevant for the current paper (see also Sects. [2] and 
[6]). It has been shown that under certain circumstances single Daphnia circle horizontally around 
a vertical artificial light shaft to which they are attracted. For high Daphnia densities a swarm 
emerges, where all Daphnia circle in the same randomly chosen direction 0, 53,01. 



The physical, biological, and chemical reasons for swarming in Daphnia in particular and in prey 
animals in general are not completely understood. Biological considerations suggest that circling 
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is the least energy consuming motion for permanently moving animals to stay as a group at a 
certain local position without frequently bumping into each other. Recent studies concentrated 
on the response of individuals as well as groups to various external influences, such as available 
food, food gradients and predator threat 25, 27, 28|, 31, 32]. In particular, 0ien [13] used methods 



from plasma kinetic theory to derive macroscopic equations - so-called fluid-dynamic equations - 
for the density of Daphnicle (Daphnia-like particle) swarms as a function of food-concentration, 
food saturation of Daphnia and a threat field of predators. The advantage of this approach is, 
that these equations can be solved analytically by approximations in linear space. However, it 
is difficult to compare these results for swarming to rotating real Daphnia swarms, as we are 
interested in. 



While 0ien [37J concentrates his description on the macroscopic level (density ond velocity de- 
stributions of the swarm), we aim at understanding swarming from a "microscopic" approach. 
That means, we derive equations of motion for individual entities - so-called agents - and inves- 
tigate the collective motion by means of a multi-agent system. This allows us to understand the 
behavior at the system level from the interactions of the entities comprising the system. 

Individual-based or agent-based modeling has turned out to be a very useful tool for modeling 
biological phenomena at various levels of organization @, 0, 0, 29]. Thus, recently different 



computer architectures have been developed to simulate the collective behavior of interacting 
agents in distributed artificial intelligence (see e.g. http://www.swarm.org/). However, due to 
their rather complex simulation facilities many of these simulation tools lack the possibility to 
investigate systematically and in depth the influence of specific interactions and parameters. 
Instead of incorporating only as much detail as is necessary to produce a certain emergent 
behavior, they put in as much detail as possible, and thus reduce the chance to understand 
how emergent behavior occurs and what it depends on. Therefore, in this paper, we follow a 
different multi-agent approach, based on Brownian agents [13] (see Sect. [3]), that - in addition 
to its computational suitability - can be also investigated by means of analytical methods from 
statistical physics and mathematics. 

The objective of the present study is to investigate the requirements on the microscopic level that 
lead to the formation of a vortex swarm on the macroscopic level. A swarm is called vortex swarm 
if animals cycle around an imaginary axis in the same rotational direction. On the global level 
(i.e. in three dimensions) one observes the emergence of a cone. As the main (cycling) motion 
takes place in the horizontal plane, we have restricted both the experimental observations and 
our simulation to two dimensions for simplicity. To find out more about physical reasons for 
vortex swarming, we first summarize some experimental observations previously reported, and 
then set-up a "minimalistic" multi-agent model to test some biologically relevant assumptions 
that may lead to the observed swarming behavior. 



2 Experimental Observations on Daphnia Motion 

In this section, we summarize some biological facts about Daphnia animals and their collective 
movement which have been reported in the literature. This information shall be used to motivate 
our swarming model in the following sections. 
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The water flea or Daphnia (see Fig.[fl) is a member of the Crustacea and are found in most fresh 
water ponds. Their body is enclosed within a carapace and their length is about 1mm to 3mm. 
Daphnia swims with a jerky motion through the water as the powerful 2nd antennae are thrust 
downward. 




Figure 1: Daphnia, Courtesy of Stephen Durr. 



Ecologically, Daphnia is extremely important in the food chains of ponds and lakes. So a sys- 
tematic investigation of their individual and collective behavior is of great interest. Ordemann et 
al. [39|, |40j, |4l|] have experimentally investigated the motion of Daphnia close to a vertical light 
shaft in both low density and high density Daphnia swarms, as we summerize in the following. 



Light source 



Water tank 




Figure 2: Sketch of the Daphnia experiments carried out by Ordemann et al. 



39 



40 



Using the experimental outline shown in Fig. [2j it was found that a single Daphnia is attracted 
to the light and starts to cycle (i.e. rotate) around the artificial light source, keeping its cycling 
direction for quite a while. In repeated experiments, however, the cycling direction may change 
to the opposite, which leads to the conclusion that single Daphnia, while rotating around the 
light beam, do not have a preferred direction of motion. 



The photograph is reprinted with the permission of Stephen Durr, 
http://www.btinternet.com/~stephen.durr/photographthree.html. For more detailed in- 

formation about Daphnia, see e.g. http://www.lander.edu/rsfox/310DaphniaLab.html, 
http: //ebiomedia. com/gall/ class ics/Daphnia/ 
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The same behavior has been observed for Daphnia swarms of sufficiently low animal density, 
where individual animals close to the light shaft cycle in both directions around the light shaft, 
frequently changing their cycling direction. 

Interestingly, the situation changes if instead of single or few Daphnia a large number of animals 
is put in the water tank. In this case, the Daphnia start again with their cycling motion, but then 
all tend to move into the same direction of motion. Prom a physical perspective, a symmetry 
break is observed, i.e., the symmetry between the two possible cycling directions (left, right 
rotation) is clearly broken toward one of the possibilities (left or right rotation). Both of these 
possibilities have the same chance to occur, but only one of them is eventually realized. The 
vortex formation as well as the symmetry break in the cycling direction are clearly self-organized 
phenomena that result from the collective interaction of many animals. In order to understand 
this in more detail, we derive a multi-agent model in the following. 



3 Agent model of biological motion in an environmental potential 

Our modeling approach is based on active Brownian particles or Brownian agents, respectively 
471 ] - Each of these agents is described by three state variables: spatial position rj, velocity Vi 
and internal energy depot ej. The first two state variables describe the movement of the agent 
and can be observed from the outside. The agent's energy depot, however, is an internal variable 
that considers the take-up of energy from the environment, the storage of energy and conversion 
of stored energy into energy of motion [48|. Provided a supercritical supply of energy from the 
environment, the Brownian agent is capable of active movement, e.g. in a preferred direction. 
The term "Brownian" refers to the fact that the agents may still be subject to fluctuations that 
are described by a stochastic force. 

The model of Brownian agents was widely discussed in different publications 11, H, 15, 48, K 



Therefore, only the basic dynamics are summarized here. For the external variables r% and 
we find the equations of motion in the form of a generalized Langevin equation: 

d d , 2 \ 



d p = v i ; J t v i = -7K) Vi - VU(r)\ ri + V2D£ i (t) (1) 

Here, for the mass m = 1 is used. Causes for the change of the variables are summarized on the 
right-hand side of the equations. The change of the agent's position, r$ is caused by the movement 
of the agent, described by the velocity Vi, that in turn can be changed by three different forces, 



explained in the following. The first term, r y(v i ), is a non-linear friction function 111. Il5j|: 



7K 2 ) =70-^2 e;(i) = 70 - - ( 2 ) 

which considers the active motion of the agent. 70 is the friction coefficient known from passive 
Brownian motion, whereas the other terms describe the influence of the internal energy depot 
ei(t), which mainly compensates this friction. Assuming that the internal energy depot relaxes 
very fast into a quasi-stationary equilibrium (adiabatic approximation), we derived an expression 
for the quasi-stationary energy depot dependent on the characteristic parameters describing its 
dynamics: qo is the influx of energy into the internal depot, which is assumed as constant here, 
c describes the loss of energy due to internal dissipation (metabolism), whereas di describes the 
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conversion rate of internal energy into kinetic energy. The nonlinear friction function has a zero 
for 

70 «2 

Active motion, i.e. \vq\ > becomes possible only for a certain supercritical take-up of energy 
from the environment, qo > 070/^2- 

The second term in eq. JT]), VC/(r), is used to describe the influences of the environment. The 
actual motion of the agent is a compromise between its active motion - which eventually would 
lead it everywhere, as long as internal energy is provided - and the environmental conditions 
which set some restrictions on this motion. The experiments described above have used a vertical 
beam of light that causes an attractive force on the Daphnia, which tend to cycle around it. In 
order to cope with this, we may choose the very simple assumption of an environmental potential 
of the form 

U(r) = V (4) 

which generates an attractive force F = — Vt7(r) = —ar towards the center, r = 0. 

Eventually, the last term in eq. ((TJ) is a stochastic force £ (assumed to be Gaussian white noise) 
of strength D, which describes the influence of random events on the agent's motion. 

The equation of motion for the Brownian agents, eq. (JTJ) , is formulated by using two dynamical 
variables, T"i(t) 3 Vi(t), as originally proposed by Langevin for the motion of Brownian particles. 
In this description, fluctuations in the environment are summarized in a stochastic force that 
changes the acceleration of the particles (according to Newton's law of motion). In the so-called 
overdamped case, one can derive from this equation the overdamped Langevin equation by as- 
suming a quasistationary velocity, v ~ 0. This results in only one equation for f, where the 
stochastic term appears with a different prefactor, y/2D r , D r being the spatial diffusion coeffi- 
cient. While such an approximation is convenient for further theoretical investigations, in the 
following we use eq. ([I]) for our computer simulations. 

Fig. [3] shows computer simulations for the active movement of a single agent, bound by an envi- 
ronmental potential, eq. , as described by eq. JT]) The result clearly indicates the cyclic motion 
round the center, which has been also observed in single Daphnia motion, as explained above. 
Running the computer simulations for single agents with different initial conditions eventually 
results in the same kind of cyclic motion, but with different rotational directions, i.e. left-handed 
or right-handed rotations. Due to stochastic influences, also changes of the direction of motion 
become possible. Thus, we may conclude that our model of Brownian agents sufficiently describes 
the observed behavior of single Daphnia. 

We now turn to the case of many, i.e. i = 1, ...,N Brownian agents, which is of importance for 
swarming. The dynamics of the multi-agent system is then described by 2N coupled (stochastic) 
equations of the form (Q]). In this case, the computer simulations shown again the characteristic 
rotational motion where, however, about half of the agents rotate clockwise, while the other half 
rotates counterclockwise. The two different cyclic directions can be clearly observed when looking 
at the angular momentum distribution, p(L), where L (for m = 1) is defined as L = r x v. As 
Fig. H] shows, this is a bimodal distribution of about equal height, indicating the both left- and 
righthanded rotational directions with the same probability. This is not surprising as long as 
independent random processes with a certain symmetry are considered. But we use this graph 
here for comparison with the results of the elaborated model, presented in the next section. 
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Figure 3: Trajectory of a single Brownian agent moving in an environmental potential, eq. (J4j, 
after t = 200. Parameters: 7 = 5.0, d 2 = 1.0, q = 10.0, c = 1.0, D = 0.005, a = 0.5, 
i.e. supercritical take-up of energy, qo > cyo/d,2- Initial conditions: {x(0),y(0)} = {0,0}, 
K(0),^(0)} = {0,0}, e(0) =0. 
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Figure 4: Angular momentum distribution p(L) of N = 1000 Brownian agents after t = 150. The 
positive or negative sign of L indicates the right- or lefthanded rotation. Parameters: qo = 10.0, 
c = 1.0 7 = 20.0, d 2 = 10.0, D = 0.001, a = 1.0. 



The simulation result also does not quite agree with the observation of high density swarms of 
Daphnia, which apparently cycle into one, i.e. the same direction. The reasons for this mismatch 
are quite obvious: in our model, we have so far only considered "point-like" agents without any 
kind of mutual interaction, whereas in real biological systems the coherent motion of the swarm 
is certainly based on local interactions between the entities. The results of the present model can 
be compared to the case of Daphnia at low density, where cycling around the light shaft in both 
directions is observed |!l| (compare Fig. [3] with Fig. 3 of Ordemann et al. and Fig. [4] with 
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Fig. la of Ordemann et al. [4lj). Thus, the question arises, which kind of interaction may lead to 
the break in the rotational symmetry, as observed in the Daphnia experiments for high animal 
density. 



4 Modeling swarming with avoidance behavior 



4.1 Interaction 



So far, different forms of global or local interactions have been introduced into swarming models. 
We mention 



3. 



1. local interactions via a self-consistent field that has been created bv the agents and in 
turn influences their further movement and/or "behavior" [1, [l3, 2^, 50, 51] - chemotactic 
response is a prominent example here. 

2. local interactions based on the coupling of the agent's individual velocity to a local average 
velocity A, @, E3, Q3, E3, 0, 3 . 



global interactions, such as the coupling of the agent's individual orientation (i.e. direction 
of motion) to the mean orientation of the swarm 0,0], or the couplin g o f the agent's 
individual position to the mean position (center of mass) of the swarm 

0,13,13 

, further 

couplings via the mean momentum or mean angular momentum or a combined set of 
invariants of motion 0, S3]. 



4. interactions based on hydrodynamic coupling between agents [IJ] 



Despite the fact that some of these models simulate coherent swarm behavior or even rotation 
of the swarm in the same direction, there is evidence that the underlying assumptions especially 
for global interactions can hardly be satisfied by biological observations, thus their biological rel- 
evance is rather questionable. Therefore, in the following section, we introduce local interactions 
between the agents that indeed match with biological reality. In particular, we focus on a special 
form of repulsive force between agents, which models avoidance maneuvers between the agents. 

Experiments on Daphnia swarming (Sect. [2]) have shown that these animals tend to cycle into 
the same direction for high Daphnia densities. We argue that the reason for this may be that 
animals try to avoid as much as possible collisions with other animals - which would occur much 
more frequently if different animals cycled into opposite directions at the same time. Thus, a 
biologically satisfactory assumption is to include avoidance behavior in our model of swarming, in 
order to test whether this would lead to the observed break in the rotational symmetry described 
above. 

Daphnia are able to sense their environment to a certain degree using their sensitive mechanore- 
ceptors fl~9| and their vision, i.e. they can detect animals approaching them from the fore, and 
then try to avoid collisions. In our models, we account for this by assuming that there is a 
short-ranged repulsive force between agents, to prevent their collisions. A similar idea was used 
by Couzin and Franks [1] to describe collective motion in ants but, different from our approach, 
they assume a hard-core repulsion in a fixed area around each agent. 
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In our model, the repulsive force results from an interaction potential V(r$) around each agent 
i that depends on its actual position, rj, and implicitly on both its actual velocity and the 
velocity of the approaching agent. In the following, two different ansatzes for such an avoidance 
potential shall be introduced. The first one discussed in the next section, has the advantage of 
mathematical simplicity, while the second one discussed in Sect. 14.31 gives smoother trajectories 
of the agents. Both ansatzes, however, lead to the same dynamic behavior of the swarm and 
therefore can be used equivalently. 



4.2 Simple Avoidance Model 

Our first approach is based on the assumption that the repulsion between two agents depends 
inversely on the Euclidian distance = \\rij\\ = \\ri — rj\\ between two agents i and j in 
two-dimensional space: 

V 1 (r ij )=-^ n G IR+, (5) 
where c is some constant. The force between two agents i and j can be calculated as 

C Tl 

fij = ~ W l (%) = (r |. )2 +n r « ( 6 ) 

Here, we have added a small offset e to the denominator, to avoid unwanted singularities if 
nj -» 0: 

4 = y/e + rij ■ r~j (7) 
In a JV agent system the total force on agent i is simply the sum over all 2 agent forces 

Fi = J2fv (8) 

The consideration of the avoidance behavior leads to a modified equation of motion, i.e. eq. (pQ) 
now reads 

jri = »i ; —Vi = - 7 K 2 ) v t -an + J2 f , L] + y/2D£(t) (9) 

Again, we have assumed a linear superposition of all these forces, which seems to us the most 
simple assuption to start with. Other assumptions for multiple interactions are of course possible, 
but hardly motivated at the current stage. 

Simulations of eq. Q with defined by eq. ((H) show a swarming behavior with an angular 
momentum distribution as in Fig. [H That means, we still find left- and righthanded rotation 
at the same time and no symmetry break. This is due to the fact that in eq. (JU) the avoidance 
behavior only depends on the distance between the two agents, i.e. the repulsion is the same 
to the front and the rear for equal distances. This assumption, however, can hardly be satisfied 
for Daphnia because they can mainly detect animals in front of them using their eye and the 
mechanoreceptors at their swimming antennas. That means that the direction of motion given 
by the velocity Vi is crucial. 

To account for this, we extend eq. ([6]) by multiplying our force with an asymmetry factor (jy that 
depends on the positions r», rj and the velocities Vi, Vj of the two agents i and j as explained 
in the following. 
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agent i 




agent j 



Figure 5: Visualization of the vectorial quantities associated with the agents. 

The prefactor has to reflect two circumstances: (i) it must increase with increasing relative 
velocity 

V Ie \ = Vij =Vi — Vj (10) 

(see Fig. [5]), because the agents would reach the point of presumable collision faster and therefore 
the force to avoid this has to be stronger; (ii) it has to consider whether two agents detect each 
other or not. This is determined by the scalar product of the relative velocity, v ie \, and the unit 
vector fy = jp^ 
each other only if 



pointing from agent i towards agent j (see Fig. [5]). Two agents i and j detect 



v Te \ ■ hj > f] > 



fill 



where r\ accounts for the angle of perception, i.e. rj = would mean a visual angle of 180° and 
77 > corresponds to a smaller one. To avoid singularities in r^, we eventually replace ||ry|| by 
eq. 0. Then, the prefactor Wy reads in its final form: 



UJ: 



V 



if (V le \ ■ Tij) >7] , T) > 

else 



(12) 



where 



(13) 



Considering the prefactor, the avoidance term eq. (JSJ) reads now 



cn 



1.1 



e \2+n ' V 



(14) 



The equations of motion are still given by eq. 
eqs. ©, (HH, (P2J will be shown in Sect. El 



Computer simulations of this extended model 
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4.3 Advanced avoidance model 



Our second approach is motivated by the observation that computer simulations of the previous 
model, while showing the correct dynamic behavior, have the visual disadvantage of abrupt 
turning maneuvers of the agents. In order to improve the visual appearance, we adopt an ansatz 
for the avoidance potential that has been originally used to model the movement of pedestrians 

V(Ri) = p-ex.pl J — J —\ (15) 

p denotes the strength and a the range of the potential, the latter being a measure of the range of 
detection. Ri is a specific function of the distance between agents, as explained in the following. 
Since all agents are moving, agent i needs to account for the space that will be occupied by all 
other agents j in the vicinity during the next time step. This space needed, depends both on the 
agent's positions rj and their velocity of motion, Vj, so Ri is a function of these. For further 
specification, we introduce the unit vector in the direction of motion of agent i, = Vi/||t>i||; n® 
is defined similarly. This allows to define a new velocity- dependent coordinate system for agent 
i, namely Pi and Xi defined by: 

Hi = Ti o x ~ou 5 X i^-Vi and (xi,Xi) = 1 . (16) 

If 5 > 0, the direction of motion of agent j is also taken into account for agent i. The a;, can 
be constructed by the orthonomalization algorithm by Gram-Schmidt. Using this coordinate 
system, the dependence of Ri on the position and velocity of agent j is now given as 



Ri = y(ri-r j ,x i y + p(r i -r j ,y i y (17) 
with a velocity- dependent function: 

/? : {ri-r^Vi) >0 

0=< t I \/n ( 18 ) 



In order to understand the meaning of the parameters 0' and A we note that, if (n — rj,y) > 
then agent i is moving away from agent j . This means that increasing the value of also means 
increasing Ri. Therefore, for both cases mentioned in eq. (fT8l) . (3' > 1 will lead to a reduction in 
the repelling force between agents i and j. On the other hand, an increase in A • vi means that 
the repelling force between agents i and j will increase if agent i is moving towards j. We note 
that these assumptions lead to an asymmetric repulsive potential V(n) around each agent. The 
potential defined by eq. (fT5j) with eq. (fT7|) can be seen in Fig. [6] for different parameters. 

Eventually, with the known repulsive interaction potential V(rj), the force between any two 
agents i and j is given as: 

= -VV(Ri) = exp ( r< - r,) (19) 

The dynamics of this model is also described by eqs. ((9j) , (EJ) . 
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Figure 6: Equipotential lines of the repulsive potential V(ri), eq. (fl~5l) for different parameters 
5 and A. Left: 6 = 0, A = 1.8, Middle: 6 = 0, A = 3, Right: 5 = 0.4, A = 1.8. The black arrow 
(in the center) indicates the agent in the origin, having a velocity of v = {0,1}. The two red 
(gray) arrows represent other agents with the same absolute value of the velocity (1) and point 
towards the origin. In the left, middle and right part of the figure, the equipotential lines shown 
correspond to the same values for the potential. 

5 Results of computer simulations of both models 
5.1 Simulation of vortex formation 

The computer simulations of both, the simple avoidance model (Sect. l4~2j) as well as the advanced 
avoidance model (Sect. l4~3l) show indeed the expected symmetry break for the swarming behavior 
in agreement with the biological observations. Spatial snapshots of a computer simulation of the 
multi-agent system with respect to avoidance behavior, together with the respective distribution 
of the angular momenta p{L) are shown in Fig. [7l The results of both models can be concluded 
as follows: 

1. On the spatial level, we observe the emergence of a coherent motion of the multi-agent 
swarm out of a random initial distribution. This collective motion is characterized by a 
unique cycling direction (either left- or righthanded rotation). 

2. We further observe the formation of a vortex which is rather similar to the Daphnia swarm 
cycling round the light beam. 

3. While in one simulation all agents cycle in the same direction, we note that in different 
simulations the cycling direction can be also opposite, i.e. there is no preferred cycling 
direction for the swarm, which also agrees with the observations of the Daphnia swarm. 
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The computer simulations for the simple avoidance model, while showing the correct dynamic 
behavior, have the visual disadvantage of abrupt turning maneuvers of the agents. To this end, 
the advanced avoidance model was introduced in Sect, 14.31 Prom the computer simulations of the 
advanced avoidance model the following improvements can be seen: 

1. The movements of the agents look smoother, more Daphnia like. 

2. Although cycling in the same direction, there are some agents that come near the center. 

3. For certain parameters a spontaneous change in the rotating direction can be observed. 
This occurs in particular if agent i takes strongly into account (5 « 0.5) the movement of 
agent j. 

In order to demonstrate the influence of the avoidance interaction on the collective motion of 
the swarm, we have conducted a computer simulation where the interaction between agents is 
"switched on" at time t = 150. Thus, in the beginning, the swarm consists of non-interacting 
agents as described in Sect. El In Fig. [8] the evolution of the angular momentum distribution 
in time is shown. In the very beginning, we find a broad distribution of p(L) centered around 
L = 0. This distribution evolves towards a clear bimodal distribution as also shown in Fig. 
IU indicating the complete symmetry between lefthanded and righthanded rotational direction. 
When the interaction potential becomes effective at t = 150, the agents start to avoid collisions 
and thus tend to move into the same direction. This can be clearly seen in Fig. [8] where p(L) 
transforms from a bimodal into an unimodal distribution after t > 150. The transformation 
period (At 30) is characterized by large fluctuations that sometimes even give the less frequent 
rotational direction a chance to take over. 

The symmetry break towards one cycling direction can be also interpreted as a process of con- 
sensus formation in groups of animals as discussed by Couzin et al. [6||. They have developed a 
model to investigate the effective leadership and decision making in animal groups in detail. To 
this end, it was assumed that a small portion of the animals is informed and therefore has a 
desired direction of motion. The model is based on a mechanism with different tunable parame- 
ters, namely the number of informed group members and the weight they give to their preferred 
direction. Through extensive computer simulations the authors find that (a) for a given group 
size the accuracy of group motion increases asymptotically as the proportion of the informed 
individuals increases, and (b) the larger the group size, the smaller the proportion of informed 
individuals needed. Couzin et al. [(| also investigate how individuals achieve consensus about the 
majority direction in case of different informed groups with competing directions. A feedback 
mechanism about adjusting the weights is proposed, to resolve this conflict. 

We like to point out that our model does not need such tunable model parameters to achieve 
consensus about the moving direction. In fact, each individual in our model, while moving, has 
only the desire to avoid collisions with other animals. Thus, the collective motion emerges as 
a self-organized phenomenon. To show that the transition from bimodal to unimodal motion is 
inherent in our model, in the following section we show how this transition occurs with increasing 
swarm size. 
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Figure 7: Spatial snapshots (left) and distribution of angular momentum p(L) (right) for 
a multi-agent system (N = 20) at three different times: (a) t = (b) t = 8 and (c) 
t = 55. The length of the arrows indicates the velocities. Initial conditions: {xj(0); y,(0)} G 
[-1.5; +1.5], {vf(0),v^(0)} = {0,0}, ^(0) = 0, parameters: 7 = 20.0, d 2 = 20.0, q = 10.0, 
c = l.o, D = 0.005, a = 1.0 p = 0.8, a = 0.1, 5 = 0.0, A = 10.0, 0' = 2. A video of the computer 
simulations can be viewed at http://intern.sg.ethz.ch/publications/2005/web-ms.html. 



5.2 Swarm size dependence of vortex formation 

So far, we have shown that our model can in fact reproduce the observed cycling behavior of 
Daphnia swarms. This means our proposed avoidance behavior on the microscopic scale indeed 
leads to the symmetry break on the macroscopic scale, whereas for single agents a symmetric 
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Figure 8: Density plot (indicated by the gray scale) of the angular distribution p(L) versus time. 
The interaction between agents is turned on at t = 150 (compare to Fig. [4]). For the parameters 
and the setup see Fig. [7l 

distribution of the angular momentum is found, in agreement with the biological observation of 
Daphnia. 

This leads us to the question whether there is a critical swarm size at which the emergence of 
a vortex can be observed. In order to investigate this, we have conducted extensive computer 
simulations of our model with a fixed set of parameters, but different swarm sizes. We point 
out that the realization of a vortex swarm strongly depends on the parameters of the model, 
in particular how much the symmetry break is enforced by the avoidance potential. Thus, any 
conclusion about a critical swarm size for vortex formation, drawn in the following, is valid only 
for the particular parameter setting. This holds also for the scaling relation discussed below. 

Secondly, we note the strong dependence of the vortex formation on stochastic influences. I.e., 
whether or not a vortex is formed, how long it takes for the establishment of a common cycling 
direction, and what this direction will be, is affected by stochastic fluctuations, which play a 
considerable role especially for small numbers of agents. 

To compensate for this, we have measured the angular momentum distribution of the swarm 
only after a sufficient time, t = 300, where a common cycling direction was established in all 
cases. This, however, does not mean that all agents follow the same direction at that particular 
time. Further, there could still be large fluctuations afterwards (as can be clearly observed in 
the computer simulation video mentioned in Fig. [7}. To account for this, we have monitored 
the angular momentum distribution over the next 50 time units, i.e. between t = 300 — 350, 
and have averaged over that time interval. We note that t = 300 does not mean simulation 
steps, but physical time, where the simulation interval was chosen as At = 5 • 10 -4 , i.e. t = 300 
corresponds to 6 • 10 5 simulation steps and the distribution was averaged over the next 10 5 
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simulation steps. It is obvious that due to the pairwise interaction of the agents by means of 
the avoidance potential, the computational effort for each simulation step also increases with 
iV 2 . Further, we have averaged the results of computer simulations over 80 runs, from which we 
calculated the sample standard deviation. 

Fig. M shows a function of the mean fraction of agents with a particular angular momentum 
dependent on the swarm size N, 

- * m 

F L+ (N) = \x L+ -x L _\ = \2x L+ -l\ ; XL+ = -—^Y,J2 N L+( t + nAt ) (2°) 

k=l n=0 

where N^)(t) is the number of agents found with a positive angular momentum in simulation 
(k) at time t, For the simulations s = 10, t = 300, At = 5 • 10~ 4 and m = 10 5 were chosen. 




Figure 9: Mean fraction of agents with positive angular momentum Tl+{N) 7 eq. (j20j) . vs. swarm 
size N . The results are averaged over 80 simulations for each data point. The fitted curves are 
given by Eqs. (|2T1) . (|22l) (black) and Eqs. ([23]) . ([24]) (red: saturation curve). Parameters: qo = 10.0, 
c = l.o 7o = 20.0, d 2 = 10.0, D = 0.001, a = 1.0, p = 0.1, a = 0.1 , 5 a = 0.01, A = 10.0. 

If a swarm has a bimodal angular momentum distribution as shown in Fig. IH xl+ would be 0.5 
and J r L_|_(A r ) = 0, i.e. no common cycling direction has been established, that is followed by a 
majority of agents. On the other hand, a clear unimodal distribution as e.g. observed in Fig. [7j 
would lead to either xl + = or 1, i.e. Fl+^N) = 1. 

As Fig. [9] indicates, the emergence of a common cycling direction occurs between swarm sizes 
of 4 to 18 agents (for the given set of parameters). Even for large swarm sizes, the common 
cycling direction is not followed by all agents, a noticeable minority fraction still cycles its own 
way. This, however, should be not considered as a drawback of the model, in fact it makes it 
much more realistic, as the computer simulation video also shows. This behavior also agrees 



15/124] 



Robert Mach, Frank Schweitzer: Modeling Vortex Swarming In Daphnia, 
Bulletin of Mathematical Biology 69 (2007) 539-562 



with observations of Daphnia swarms at high density. In particular, animals at the border of the 
swarm still do not follow the cycling direction of the majority. 

We note that the occurence of a common cycling direction is not an abrupt transition, but 
a gradual one. One can argue that this resembles a kinetic phase transition analogous to the 



continuous phase transition in equilibrium systems [541 ] . i.e. 



F L+ (N) =co(N-N c ) K (21) 

where N c is the onset of the transition and k is the critical exponent. For their model Vicsek 
et al. [3] have determined k as 0.35. We have tested this finding against our simulation results, 
as shown in Fig. [9] and found for the least square fit and standard errors in the region N G [4, 20]: 

N c = 3.240 ± 0.328; c = 0.422 ± 0.029; k = 0.270 ± 0.025 (22) 



On the first glimpse, there is a close similarity between Vicsek et al. [54| and our findings with 
respect to the value of the critical exponent k. I.e., we emphasize that vortex formation under 
repulsion also occurs at a critical point. 

However, as Fig. [9] clearly shows, the scaling suggested by Vicsek et al. [sij only holds for a 
restricted range of N - which by the way was already the case in the original paper (cf. Fig 2b 
of Vicsek et al. jE^]). A much better fit of the observed transition is provided by the following 
function: 

N-Nc~ 



(23) 



1 + exp 

with the best-fit parameters and standard errors: 

N c = 4.496 ± 0.136; r = 3.651 ± 0.161; « = 0.880 ± 0.004 (24) 

As shown in Fig. [9j this function fits the simulated data extremely well. The only outlier for 
iV = 3 is due to the fact, that for 3 agents at least 2 agents go in the same direction. That means 
that the probability of breaking the symmetry is a lot higher in this particular case. 



In addition to the depencence on the swarm size or density, Vicsek et al. [54| also discussed a 
phase transition with respect to noise, rj, in the alignment of the individuals and suggested a 
scaling similar to Eqn. (|21~T) . As the critical value, rj c , depends on the lattice size of their computer 
simulations, there is no direct comparison with this result possible. However, we note the obvious 
similarities to other types of structure formation which also occur only below a certain critical 
temperature. This holds also for our model where a common cycling direction can emerge only 
below a critical noise level. As a difference to the investigations of Vicsek et al. [54], noise in 
our model enters the equation of motion, Eqn. ([9]), whereas the avoidance interaction is without 
noise. 



6 Model testing of avoidance maneuvers 

So far, reasonable assumptions about local interactions of Daphnia (such as local repulsion) have 
been taken into account. But it still remains to be tested empirically and theoretically whether 
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avoidance maneuvers of Daphnia are really responsible for the symmetry break observed in high 
density swarms. 

To gather evidence on the microscopic dynamics of our model, we used data provided by Anke 
Ordemann [H] about avoidance maneuvers of pairs of Daphnia approaching each other in a 
horizontal plane. Snapshots of a typical avoidance time series are shown in Fig. [TUJ A thorough 
analysis of similar experiments would allow to determine the parameters in our model, but the 
time series presented was the only data made available to us. Nevertheless, this sequence already 
enables us to qualitatively compare the experimental findings with our theoretical model of 
avoidance behavior. 
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Figure 10: Experimental observation of an inter-animal avoidance maneuver of two Daphnia in 
the horizontal plane In our representation of the available data, each arrow represents the 
spatial orientation (head and tail) of the Daphnia at successive times (every 0.2 sec). In the 
sequence, same arrow colors (gray scales) correspond to the same time. 

To this end, we have simulated the avoidance maneuver of two agents in a situation similar to 
the one shown in Fig. [TO) From the simulation results shown in the left part of Fig. [XT] it can be 
clearly seen that the two agents avoid each other not in a symmetrical way, but rather similar 
to the experimental findings. One may argue that this is due to the additional influence of the 
environmental potential which may more affect the agent closer to the origin (0,0). In order to 
test this, we have simulated the same situation without the avoidance interaction potential, but 
just the environmental potential. The result, shown in the right part of Fig. [TTJ clearly indicates 
that the environmental potential has only very little affects the agent's trajectories and thus 
cannot be responsible for the realistic avoidance behavior shown in the left part of Fig. [TTJ 

Comparing the simulation of the model with the above experimental observations in Daphnia, it 
becomes obvious that inter-animal avoidance maneuvers similar to the ones incorporated in the 
current Brownian agent model are most likely involved in the mechanism that causes symmetry 
breaking for high Daphnia density and leads to the observed vortex swarms. 
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Figure 11: Left: Simulation of an inter-agent avoidance maneuver in the horizontal plane. Each 
arrow represents the spatial orientation (drawn from the velocity vector) of the agent at successive 
times (every 0.2 time units). In the sequence, same arrow colors (gray scales) correspond to the 
same time. Parameters: q = 10.0, c = 1.0 70 = 20.0, d 2 = 10.0, D = 0.001, p = 0.2, a = 0.2 
, 5 a = 0.1, A = 10.0, a = 0.1, Right: Simulation with the same setting but without avoidance 
interaction, p = 0.0. 

7 Discussion 

In this paper, we tried to understand the vortex swarming behavior observed in Daphnia from 
rather minimalistic assumptions that, however, should have a clear biological relevance. To this 
end, we introduced a multi-agent model based on the concept of Brownian agents. Different 
from other modeling approaches which are based e.g. on cellular automata, our model considers 
both continuous time and space. Further, in addition to external and interaction forces effecting 
the agent's behavior, we consider stochastic influences resulting e.g. from random events or 
fluctuations in the environment. As a difference to multi-agent models developed in artificial life, 
our model further allows for a more specific discussion of the different parameters affecting the 
agent's motion. It can be also seen as an alternative to some of the self-propelled particle models 
17, III, 54] where alignment between individuals is explicit. In our model, we demonstrate that 



vortex swarming does not require an explicit alignment rule for individuals but can be achieved 
through a combination of attraction and repulsion. This insight allows us to reduce assumptions 
about the cognitive requirements of swarming organisms. 

At the end, we give a critical discussion of the Brownian agent model in the biological context of 
Daphnia swarming, in order to point out its limitations and to give some hints for experimental 
verifications. 

The equation of motion for the Brownian agents is a generalized Langevin eq. ([T]), which models 
the motion of the animals as continuous in space and time. The movement of single Daphnia, 
however, consists of a series of jumps as can be also observed in experimental observations. 
This can at least be partially covered by using a larger time step in our simulation, resulting in 
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larger jumps during each time interval. We argue that this would not change the overall dynamic 
behavior of our model. Thus, we kept the continuous approximation of the discrete jumps and 
the small time step to avoid numerical instabilities. 

Our theoretical description is based on the assumption of a quasistationary energy depot of the 
agents. This implies a constant take-up of energy, qo, whereas Daphina feed in cycles while they 
do not move actively Thus, a switch-like change in behavior seems to be more appropriate. This 
can be covered in our model by assuming the take-up rate of energy a time-dependent function, 
q(t), and replacing the continuous acceleration term in eq. (pQ), (fee^i), by a more complex term 
that reflects a switch-like change between starvation and active motion. 

Response to external forces and force-based interaction between the agents are two basic ingre- 
dients of our model of Brownian agents. Interaction of biological entities, on the other hand, is 
often driven by internal and behavioral reasons. This certainly sets the limits for our approach 
which is inspired by, and based of, physical considerations using minimal assumptions. Physics, 
per se, has no concept of "behavior" based on internalized motivations of an agent. Thus, in 
our model every kind of internal "driving forces" has to be externalized by assuming that the 
agent behaves as it would follow an external force that leads it to the desired behavior. Such an 
approach does not claim that these forces really exist in the outer world, it is rather a convenient 
modeling formalism that allows to apply the concepts of physics to the much more complex be- 
havior of animals. An alternative could be concepts from artificial intelligents, where the internal 
dynamics of agents is modeled explicitely (the so-called BDI (belief-desire-intention) agents, for 
example, can have their own internal world view). 

In our model of Brownian agents, we could demonstrate the emergence of a vortex swarm based 
on local asymmetric interactions of the agents. The global dynamic behavior was also found for 
other models (see also Sect. I4.1| ) exploiting other mechanisms of symmetry breaking. We do not 
want to argue here about the most simplest one - at the end, one has to find a compromise 
between "simple" and biologically "satisfied". Many of the proposed mechanisms - such as the 
hydrodynamic coupling - still lack an experimental justification. Our avoidance model, on the 
other hand, could be at least visually tested by comparing the experimental observations of two 
Daphnia approaching each other with a similar situation from the computer simulations. 

Noteworthy, in our model the emergence of a vortex swarm is not enforced by the alignment 
of the agents, as used in other models. Instead, we have included only the simple assumption 
of collision avoidance. The specific form of our avoidance potential penalizes mainly the head- 
head collision and thus promotes the dispersion of the agents (and not the alignment). As a 
consequence, to obtain a vortex swarm, an additional attraction force is needed, which in the 
considered case results from the environmental potential (attraction towards the center). This is 
not a drawback of the model, but justified by the real experimental observation. As shown by 
Schweitzer et al. [Hi], under certain circumstances the effect of an external parabolic potential is 
also equivalent to an attractive force between the agents. The emergence of a vortex swarm can 
be seen as a dynamic compromise between three different requirements: active motion (to keep 
the agents moving), asymmetric repulsive forces (e.g. to avoid collisions) and attractive forces 
(either enforced by a local agent-agent attraction, or an external potential). While the first two 
requirements alone would simply lead to a dispersion, the latter one results in a compression (or 
confinement) of the swarm. From a physics perspective the vortex swarm is a stable attractor 
of the multi-agent dynamics; from an "economics" perspective one could think of global utility 
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maximization balancing out all individual requirements (such as avoiding discomfort from both 
collisions and separation). 

A noticeable advantage of our model is that it does not break down for small swarm sizes, so 
it can be used to simulate both single animals and swarms of low and of high density, whereas 
other models mainly concentrate on the dynamics of reasonably sized swarms. This allows in 
particular to investigate the transition from the uncorrelated (bi-directional) rotation of single 
agents to the correlated vortex formation of the swarm, as shown in Fig. [9J 

So far, we have shown a qualitative agreement between our computer simulations and some 
experimental observations in Daphnia. A quantitative verification, or even a prediction of Daphnia 
behavior under different real conditions is still missing. In this respect, our model is not different 
from other theoretical models proposed. A first step towards verification involves the experimental 
determination of the parameters. In this paper, we can just propose some ideas: 

• Neglecting the jump-like motion of single Daphnia and using the approximation of a (quasi- 
stationary velocity, one should be able to estimate an average velocity of cycling Daphnia. 
This can be directly related to the stationary velocity Vo, eq. j3]), that enters the equation 
of motion used in the model. 

• The parameters determining the avoidance behavior may be estimated by a direct com- 
parison between the experimental observation, Fig. [10] and the computer simulation, Fig. 
ITTl Obviously, the sole event shown in this paper, is not sufficient for that; so, we would 
expect further experimental investigations here. 

• Another way to (indirectly) estimate the parameters of the avoidance potential is via the 
local swarm density, or the (average) spatial extension of a swarm of given size. These are 
determined by our model parameters and could be possibly compared with experiments 
(varying both the strength of the light beam and the swarm size). 

• By varying the swarm size, one can also experimentally test the onset of the vortex swarm- 
ing, and compare this to the respective computer simulations shown in Fig. [9] As mentioned, 
the transition range towards the vortex formation in the model strongly depends both on 
swarm size and on the parameters characterizing the avoidance potential. So, in addition 
to the observation of avoidance maneuvers of two animals, this yields a macroscopic veri- 
fication for the parameters of the avoidance potential. 

Eventually, we want to point out some situations where the model could make predictions about 
Daphnia behavior which may be tested experimentally. The model uses the assumption of an 
environmental potential, that in the current investigation results in an attraction toward the 
center, this way considering the influence of the vertical light beam on the Daphnia. The real 
influence could be tested by producing Daphnia mutants that are insensitive to light. Then we 
expect no attractive force, and hence no vortex swarming, in agreement with the explanations 
above. Also, one could think of two different vertical light beams in the water tank, at a certain 
distance. Starting from a homogeneous spatial distribution of agents, the model would predict 
the occurrence of two different vortex swarms around the two centers, each probably having its 
own rotational direction, as long as the distance between the light beams is large enough. For 
smaller distances (where the critical distance may be also a function of the agent density in the 
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system), we would expect from the model interferences at the boarder between the two rotating 
swarms, which would lead to additional couplings and thus maybe to a synchronization of the 
rotational directions. The situation of the two separated vertical light beams would also allow us 
to test whether the attraction of the Daphnia is limited by some maximum local density of the 
swarm, as suggested by our model. In this case, we would find in the experiments two distinct 
swarms, whereas without saturation effects, one could possibly find one swarm only, with a much 
higher density. 

A remaining question is weather the model in the current form is also applicable to other species. 
Vortex swarming, as we have pointed out, is a widely spread phenomenon observed also in fish, 
or bacteria. We argue that the principle (qualitative) features of vortex swarming are covered by 
our model, as long as (local or global) attraction and asymmetric repulsion (e.g. via avoidance 
maneuvers) play a considerable role. This may hold for fish, but probably not for gliding bacteria 
like Paenibacillus, which also show vortex formation In the latter case, adhesion forces may 
play a much more important role as one can also deduce by looking at the sharp external 
boundaries of the bacterial swarm. However, it would be still possible to adapt the model of 
Brownian agents also for this situation by including other terms of local interaction. 
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